!*********************************************************************************************************************************
!Main program unit 
!*********************************************************************************************************************************
!Main steps"
!(1) Initialize the parameter vector and parameter matrix used for log likelihood maximization
!(2) Read in data
!(3) Evaluate likelihood at perturbed parameter vector (at the estimated parameter values, also compute at this step asymptotic standard errors)
!(4) Using amoeba to do log likelihood maximization
!
!NOTE: to run standard errors:
!(i)   set comp_sterr=1 in the subroutine SUBROUTINE matrix_initial(pr) in E_mod_param_initial
!(ii)  set est_ind = 1 for each estimated param
!(iii) set the variable num_param1 in the main program unit equal to the total number of estimated params
!(iv)  check that totpop_trim equals the actual # of individuals whose histories enter the likelihood 
!(v)   computest_err = 1 here below in the main program unit

PROGRAM main


USE mod_types
USE mod_param_initial
USE mod_parameters
USE mod_grid


IMPLICIT NONE



!---------------------
!(1) Usual parameters
!---------------------
INTEGER, DIMENSION(4):: sumtemp
INTEGER, DIMENSION(4,4):: arraytemp

REAL(dp), PARAMETER:: ftol = 1.0d-10
INTEGER, DIMENSION(phi_grid_size):: task
REAL(dp):: phi_1_star, phi_2_star 

!input to amoeba (initialization)
REAL(dp), ALLOCATABLE:: param_matrix(:,:)
  
!arrays storing the observed variables
INTEGER, DIMENSION(:,:), ALLOCATABLE:: level_o, perf_o, sal_o
INTEGER, DIMENSION(:), ALLOCATABLE:: edu, age, year
INTEGER, DIMENSION(:,:), ALLOCATABLE:: ind
INTEGER, DIMENSION(:,:), ALLOCATABLE:: counter
INTEGER, DIMENSION(:), ALLOCATABLE:: indelete
INTEGER, DIMENSION(:), ALLOCATABLE:: sumind, countind
INTEGER:: totind

!input to amoeba:likelihood values for rows of param_matrix
REAL(dp), ALLOCATABLE:: F(:)   
REAL(dp), DIMENSION(CAP_N):: loglike  
INTEGER:: i101, m1, m2, m3, iter, totsum, m4		       
!CHARACTER(dp):: hour
REAL(dp):: etime, try


!-------------------------------
!(2) Standard error parameters
!-------------------------------
!sample size
INTEGER, PARAMETER:: totpop =1430
!# of inds actually entering the likelihood
INTEGER, PARAMETER:: totpop_trim = 1426
INTEGER, DIMENSION(1430):: indicator_delete
INTEGER, DIMENSION(4):: indiv_delete
INTEGER, DIMENSION(5):: indicator_ex
INTEGER, DIMENSION(5):: x_ex
INTEGER, DIMENSION(3):: x_short

!matrices for the BHHH method of estimating the inverse information matrix
REAL(dp),DIMENSION(:,:), ALLOCATABLE:: INVINFO, INVINFO2, INFO, INFO2
!matrices for estimated the gradient of the log-likelihood
REAL(dp), DIMENSION(totpop):: GHAT
REAL(dp), DIMENSION(:,:), ALLOCATABLE:: GHAT2, GHAT3
REAL(dp), DIMENSION(:,:), ALLOCATABLE::GHAT3PRM 
!declare variables for standard errors and tstat
REAL(dp), DIMENSION(:), ALLOCATABLE:: SE, TSTAT
!deviation parameter
REAL(dp), PARAMETER:: CAP_ST= 0.025_dp 
REAL(dp), DIMENSION(:), ALLOCATABLE:: theta
REAL(dp), DIMENSION(1426):: loglike_short
!# of parameters actually estiamted
!INTEGER, PARAMETER:: num_param1 = 101
INTEGER, PARAMETER:: num_param1 = 75
INTEGER:: computest_err



!------------------------------------------------------------------------------------
!REMEMBER: modules with a SAVE attribute are MOD_PARAM_INITIAL and MOD_PARAMETERS
!          they should always be included in any module which declares any variables,
!          to avoid that they become undefined (Fortran 90: page 110, page 363)
!------------------------------------------------------------------------------------
OPEN(unit=3000, file = 'estimation_output.txt', status = 'replace')
OPEN(unit=3500, file = 'standard_errors_t.txt', status = 'replace')
WRITE(3000,*) ""
WRITE(3000,*) "MLE of the Three-Task Model with Wages (N=1430 (tot:1426))"
WRITE(3000,*) "(Types=4, Biased CE)"
WRITE(3000,*) "ONLY THE FIRST 7 PERIODS (fit until t=7)"
WRITE(3000,*) ""


PRINT*, ""
PRINT*, "MLE of the Three-Task Model with Wages (N=1430 (tot:1426))"
PRINT*, "(Types=4, Biased CE)"
PRINT*, "ONLY THE FIRST 7 PERIODS (fit until t=7)"
PRINT*, "" 


!*************************************************************************************
!(1) This subroutine fills a parameter with the current time as the string  hh:mm:ss
!*************************************************************************************
!OPEN(unit = 2000, file = 'Time.txt', status = 'replace')
!CALL TIME(hour)
 CALL cpu_time(etime)
 WRITE(3000,*) "Starting Time: (",(etime),")"
 PRINT*, "Starting Time: (",(etime),")"
!CLOSE(unit = 2000)


!**************************************
!(2) Initialize the parameter vector
!**************************************
!it is contained in the subroutine mod_parameters
!CALL parameter_generation(x_vector)

!(*i*) set afresh the file in which to store the actual values of the parameters
OPEN(unit=3100, file ='track_actual_param.txt', status = 'replace')
CLOSE(3100)

OPEN(unit=3100, file ='track_actual_param.txt', position="append")

!(*ii*) it is contained in the subroutine mod_param_initial
CALL param_initialization



!*****************************************************************************************
!(3) Allocate the matrix of parameters and the corresponding vector of likelihood values
!*****************************************************************************************
ALLOCATE(param_matrix(num_param+1,num_param), F(num_param+1))


WRITE(3000,*) "# of Parameters to be estimated: ", num_param
WRITE(3000,*) ""
PRINT*,       "# of Parameters to be estimated: ", num_param



!**************************************************************************************************************
!(4) Initialize the matrix P input to the amoeba subroutine (the subroutine is contained in mod_param_initial) 
!**************************************************************************************************************
CALL matrix_initial(param_matrix)


!PRINT*, "Initial Matrix of Parameters to be Estimated"
!PRINT('(3F16.4)'),   ((param_matrix(i,j),j=1,num_param),i=1,num_param+1)



!********************
!(5) Output files
!********************
!(*i*) set afresh the file in which to store the values of the parameters
OPEN(unit=777, file = 'track_like.txt', status = 'replace')
CLOSE(unit=777)

!(*ii*) set afresh the file in which to store the normalized values of the parameters
!OPEN(unit=77, file =  'track_param.txt', status = 'replace')
!CLOSE(unit=77)



!***********************************************************************
!(6) Value of the likelihood at the num_param+1 points in param_matrix
!***********************************************************************
!----------------------
!(a) Read in the data
!----------------------
!---------------------------------------------------------------------------
!i. The program is reading, for each individual in this order (full sample):
!   level_o(ind,tenure), sal_o(ind,tenure), perf_o(ind,tenure)
!---------------------------------------------------------------------------


ALLOCATE(level_o(CAP_N,max_t), sal_o(CAP_N, max_t), perf_o(CAP_N,max_t))
ALLOCATE(edu(CAP_N), age(CAP_N), year(CAP_N))


m1 = 0

!For the original estimation: use file fortran_sample.txt
!For replication material: use file fortran_sample_artificial.txt
OPEN(unit=2, file = 'fortran_sample_artificial.txt', position='rewind') !use fortran_sample_artificial.txt
m1 = 0
DO m1 = 1, CAP_N
  READ(2,*)   level_o(m1,1),  sal_o(m1,1), perf_o(m1,1), level_o(m1,2),  sal_o(m1,2), perf_o(m1,2), &
              level_o(m1,3),  sal_o(m1,3), perf_o(m1,3), level_o(m1,4),  sal_o(m1,4), perf_o(m1,4), &
	          level_o(m1,5),  sal_o(m1,5), perf_o(m1,5), level_o(m1,6),  sal_o(m1,6), perf_o(m1,6), &
	          level_o(m1,7),  sal_o(m1,7), perf_o(m1,7), level_o(m1,8),  sal_o(m1,8), perf_o(m1,8), &
	          level_o(m1,9),  sal_o(m1,9), perf_o(m1,9), level_o(m1,10), sal_o(m1,10),perf_o(m1,10), &
			  level_o(m1,11), sal_o(m1,11),perf_o(m1,11), edu(m1), age(m1), year(m1)
  END DO 
CLOSE(unit=2)



!-------------------------------------
!ii. Keep track of original histories 
!-------------------------------------
OPEN(unit=444, file="trimmedsample.txt", status="replace")
m1 =0
DO m1 = 1, CAP_N
  

WRITE(444,'(37I6)') 1, level_o(m1,1), perf_o(m1,1), & 
           level_o(m1,2), perf_o(m1,2), level_o(m1,3), perf_o(m1,3), &
		   level_o(m1,4), perf_o(m1,4), level_o(m1,5), perf_o(m1,5), &
		   level_o(m1,6), perf_o(m1,6), level_o(m1,7), perf_o(m1,7), &
		   level_o(m1,8), perf_o(m1,8), level_o(m1,9), perf_o(m1,9), &
		   level_o(m1,10), perf_o(m1,10), level_o(m1,11), perf_o(m1,11), edu(m1), age(m1), year(m1), &
		   sal_o(m1,1), sal_o(m1,2), sal_o(m1,3), sal_o(m1,4), sal_o(m1,5), sal_o(m1,6), &
		   sal_o(m1,7), sal_o(m1,8), sal_o(m1,9), sal_o(m1,10), sal_o(m1,11)


END DO

CLOSE(unit=444)


ALLOCATE(INVINFO(num_param1,num_param1), INVINFO2(num_param1,num_param1), INFO(num_param1,num_param1), INFO2(num_param1,num_param1))
ALLOCATE(GHAT2(totpop_trim,num_param1), GHAT3(totpop_trim,num_param1))
ALLOCATE(GHAT3PRM(num_param1,totpop_trim))
ALLOCATE(SE(num_param1), TSTAT(num_param1),theta(num_param1))


!----------------------------------------------------------
!(b) Initialization of the likelihood (loop on parameters)
!----------------------------------------------------------
!OPEN(unit=2000, file = 'probability_level.txt')
 !WRITE(2000,'(1A15,4A5,1A15)') 'SUMVSIM','i1','i2','i3','i4','Prob_level'
!===========================
!loop on parameters
 DO i101 = 1, num_param+1
!===========================
  WRITE(3000,*) ""
  !CALL TIME(hour)
  CALL cpu_time(etime)
  WRITE(3000,*) "Initialization of Likelihood for Row:",i101, "(",(etime),")"
  
 !IF (i101==5) THEN
 !  PAUSE
 !  PRINT*, 'Fifth likelihood row', i101
 !END IF

  WRITE(3000,*) ""
  !CALL TIME(hour)
  CALL cpu_time(etime)
  PRINT*,       "Initialization of Likelihood for Row:",i101, "(",(etime),")"
  

  OPEN(unit = 777, file = 'track_like.txt', position="append")
  !OPEN(unit = 77, file = 'track_param.txt', position="append")
  WRITE (777,*) "Initial Value of Likelihood at Row: ", i101
  WRITE(3100,*) "Initialization of Likelihood for Row: ", i101
  WRITE(3100,*) ""
  !WRITE (777,*) " "
  !WRITE (77,*)  "Initial Parameters, Row:", i1 
  !WRITE (77,'(5F15.10)') (param_matrix(i1,1)), (param_matrix(i1,2)),&
        !(param_matrix(i1,3)), (param_matrix(i1,4)), (param_matrix(i1,5))
  !WRITE (77,'(5F15.10)') (param_matrix(i1,6)), (param_matrix(i1,7)), &
        !(param_matrix(i1,8)), (param_matrix(i1,9)), (param_matrix(i1,10))
  !WRITE (77,'(4F15.10)') (param_matrix(i1,11)), (param_matrix(i1,12)), &
        !(param_matrix(i1,13)), (param_matrix(i1,14))
  !WRITE (77,*) " " 
  CLOSE(unit=777)
  !CLOSE(unit=77)




!------------------------------------------------------------------
!(c) Evaluation of the likelihood at the initial parameter vector
!------------------------------------------------------------------
  CALL likelihood_function(param_matrix(i101,:),F(i101),loglike(:))
  

     
  OPEN(unit=777, file='track_like.txt',position="append")
  WRITE(777,*) (F(i101))
  WRITE(777,*) ""
  CLOSE(unit=777)
 



!**************************************************************************************************
!(7) Compute asymptotic std errors. The asymptotic variance of the estimates is computed using the
!    BHHH estimator, proposed by Berndt, Hall, Hall and Hausman in '74, based on the outer product 
!    of the score, the gradient of the log-likelihood at the estimated parameter vector 
!    (check smaller values of the percentage change in parameters (CAP_ST= 0.025_dp), make
!     sure this change reflects the same absolute difference in the original and perturbed
!     value of a parameter, and double check the same value is imposed on the bump in 
!     SUBROUTINE matrix_initial(pr) in E_mod_param_initial)
!***************************************************************************************************
!indicator function for computing standard errors
computest_err = 1


!PRINT*, 'num_param =', num_param
!PAUSE


!---------------------------
IF (computest_err==1) THEN
!---------------------------
PRINT*, "Start the computation of the Info Matrix"


!---------------------------------------------------------
!(a1) Delete individuals not included in the likelihood  
!---------------------------------------------------------
indicator_delete = 0

OPEN(unit=4000, file = 'ind_delete.txt', position='rewind')
m1 = 0
DO m1 = 1, 4
  READ(4000,*) indiv_delete(m1)
  !WRITE(3000,*) ind_delete(i1)
  indicator_delete(indiv_delete(m1)) = 1
END DO

CLOSE(unit=4000)

!--------------------------------------------------------------------------------------------------------
!(a2) Create a shorter likelihood vector which does not include individuals recorded in 'ind_delete.txt'
!--------------------------------------------------------------------------------------------------------
m1 = 0
m2 = 0
m3 = 0



DO m1 = 1, 1426
 
!checking the value of indicator_delete ov
IF(indicator_delete(m1 + m3) == 0) THEN
  
  loglike_short(m1) = loglike(m1+m3) 
 
 ELSE IF(indicator_delete(m1+m3) == 1) THEN

  IF (m1+ m3 + 1 <= totpop) THEN

   DO m2 = m1+m3+1, totpop
    IF (indicator_delete(m2) == 0) THEN
     loglike_short(m1) = loglike(m2)
	 m3 = m2 -m1 -m3 +m3  != no. of inds to skip = new + old = m2 - (m1+m3) + m3
	 EXIT
    END IF
   END DO
 
  ELSE IF(m1+ 1+ m3 > totpop) THEN
   EXIT
  END IF

 END IF


END DO



OPEN(unit=4004, file = 'loglike_short.txt', position='rewind')

DO m1 = 1, 1426
 WRITE(4004,*) loglike_short(m1)
END DO

CLOSE(unit=4004)



!---------------------------------------------------------------------
!(a3) Compute original and perturbed likelihood for each observation
!---------------------------------------------------------------------
PRINT*, "Compute the original likelihood for each observation"

m1 = 0
m2 = 0


WRITE(3000,*) 'GHAT (original like for trimmed pop)'
!================
IF(i101==1) THEN
!================
   !construct a matrix of likelihood values at the estimated parameter vector of all inds
   DO m2 = 1, totpop_trim
     GHAT(m2) = loglike_short(m2)
     WRITE(3000,*) GHAT(m2)
   END DO


!-------------------------------------------------------------------------------------
!(b) Likelihood computed at the perturbed parameter vector (the numerical gradient 
!    of the likelihood is calculated using the likelihood evaluated at the perturbed 
!    values of the parameters, that is, their initial value + bump. So for every row 
!    different from the first one, the likelihood is evaluated at p_vector + bump)
!-------------------------------------------------------------------------------------
PRINT*, "Compute the perturbed likelihood for each observation"


!=====
 ELSE
!=====

 WRITE(3000,*) 'GHAT2 (perturbed like for trimmed pop)'
  
  !construct a matrix of perturbed likelihoood values of all inds
  DO m2 = 1, totpop_trim
    GHAT2(m2, i101-1) = loglike_short(m2) 
	!WRITE(3000,*) GHAT2(m2, i101-1)
  END DO

  PRINT*, 'ok'

!=======
 END IF
!=======

!---------------------------------------------------------
END IF !case if for when the st errors are to be computed
!---------------------------------------------------------
!============================
 END DO  !loop on parameters
!============================



!---------------------------------------------------------------------------------------
!(c) For each observation, estimate the gradient
!    a. the vector of estimated parameters is given by the first row of the P matrix
!    b. add all other parameters  
!    c. REMEMBER to set the vector of bumps, cap_st
!----------------------------------------------------------------------------------------
!---------------------------
IF (computest_err==1) THEN
!---------------------------

PRINT*, "Create the matrix of differences"

OPEN(unit=4100, file ='theta.txt', status = 'replace')


theta(1)  = alpha_1 
theta(2)  = alpha_2 
theta(3)  = xi_1
theta(4)  = xi_2
theta(5)  = xi_3   
theta(6)  = gammapast11_0
theta(7)  = p_1
theta(8)  = gammap_01
theta(9)  = gammap_02
theta(10) = gammap_03
theta(11) = s_1
theta(12) = s_2 
theta(13) = sigma1_u
theta(14) = s1_5
theta(15) = s_3
theta(16) = s_002
theta(17) = s_003
theta(18) = sigma2_u
theta(19) = sigma3_u
theta(20) = w_75
theta(21) = w_76
theta(22) = w_77
theta(23) = w_78
theta(24) = w_79
theta(25) = gammap_04
theta(26) = d10
theta(27) = d5
theta(28) = phk2
theta(29) = d12
theta(30) = xi_14
theta(31) = p_3
theta(32) = xi_34
theta(33) =  d
theta(34) = d16
theta(35) = xi_21
theta(36) = d1_1
theta(37) = d1_2
theta(38) = cpen 
theta(39) = xi_25 
theta(40) = alpha_3
theta(41) = d20
theta(42) = d21
theta(43) = d23
theta(44) = d24
theta(45) = beta_1
theta(46) = beta_2
theta(47) = d27
theta(48) = d28
theta(49) = d30
theta(50) = d31
theta(51) = d32
theta(52) = plus
theta(53) = tencoef1
theta(54) = beta_3
theta(55) = sigma22_u
theta(56) = sigma23_u
theta(57) = sigma24_u
theta(58) = s_0032
theta(59) = s_0033
theta(60) = s_0034
theta(61) = s_0022 
theta(62) = s_0023 
theta(63) = s_0024
theta(64) = s1_52
theta(65) = s1_53
theta(66) = s1_54
theta(67) = d5t2 
theta(68) = d5t3  
theta(69) = d5t4
theta(70) = sigma12_u 
theta(71) = sigma13_u 
theta(72) = sigma14_u 
theta(73) = s3_1 
theta(74) = s3_2 
theta(75) = s3_3 


m4 = 0
DO m4 = 1, num_param1

 WRITE(4100,*) theta(m4)

END DO

CLOSE(4100)

!PRINT*, "Theta_1 is equal to: ",  alpha_1

 
!---------------------------------------------------------------------------------------------------
!(d) Compute the matrix of first derivatives of the likelihood function
!    GHAT:: likelihood values at the estimated parameter vector for all individuals in the sample
!    GHAT2: matrix of likelihood values at the perturbed parameter vector
!---------------------------------------------------------------------------------------------------
 m1 = 0
 m2 = 0
 DO m1 = 1, totpop_trim
  
  DO m2 = 1, num_param
   
   GHAT3(m1,m2) = (GHAT2(m1,m2)-GHAT(m1))/(CAP_ST*theta(m2))
    
   IF (m1 == 27 .or. m1 ==28 .or. m1==31 .or. m1==39 .or. m1==47 .or. m1==57 .or. m1==69) THEN
    PRINT*, 'ind', m1
	PRINT*, 'like', GHAT2(m1,m2)
	PRINT*, 'perturbed like', GHAT(m1)
	PRINT*, 'derivative', GHAT3(m1,m2)
   END IF

   IF( DABS(GHAT3(m1,m2)) < 1.0d-20) THEN
     GHAT3(m1,m2) = 1.0d-20
   END IF

  END DO
 END DO


   PRINT*, "For individual 1"
    m2 = 0
    DO m2 = 1, num_param
     PRINT*, 'likelihood change =', GHAT2(1,m2)-GHAT(1)
	 PRINT*, 'size of the change =', CAP_ST*theta(m2), 'in parameter', m2
     PRINT*, 'derivative of the likelihood =', GHAT3(1,m2)
    END DO

   PRINT*, "For individual 2"
     m2 = 0
    DO m2 = 1, num_param
     PRINT*, 'likelihood change =', GHAT2(2,m2)-GHAT(2)
	 PRINT*, 'size of the change =', CAP_ST*theta(m2), 'in parameter', m2
     PRINT*, 'derivative of the likelihood =', GHAT3(2,m2)
    END DO
  
   PRINT*, "For individual 3"
     m2 = 0
    DO m2 = 1, num_param
     PRINT*, 'likelihood change =', GHAT2(3,m2)-GHAT(3)
	 PRINT*, 'size of the change =', CAP_ST*theta(m2), 'in parameter', m2
     PRINT*, 'derivative of the likelihood =', GHAT3(3,m2)
    END DO

   
!-------------------------------------------------
!(e) Compute the transpose of the gradient matrix
!-------------------------------------------------
 PRINT*, "Transpose the matrix of differences"

 GHAT3PRM = TRANSPOSE(GHAT3)
 

 !checks on the original matrix of derivatives
 DO m1 = 1, num_param
  DO m2 = 1, totpop_trim
  
  IF ( DABS(GHAT3(m2,m1)) < 1.0d-30) THEN
   PRINT*, 'like diff', GHAT3(m2,m1)
   PRINT*, 'ind', m2
   PRINT*, ''
   !PAUSE
  END IF

  END DO
 END DO

 
 !checks on the transpose matrix of derivatives
 DO m1 = 1, num_param
  DO m2 = 1, totpop_trim
  
  IF ( DABS(GHAT3PRM(m1,m2)) < 1.0d-30) THEN
   PRINT*, 'transpose of the like diff', GHAT3PRM(m1,m2)
   PRINT*, 'ind', m2
   PRINT*, ''
   !PAUSE
  END IF

  END DO
 END DO



!------------------------------------------------------
!(f) Pre-multiply the gradient matrix by its transpose
!------------------------------------------------------
 PRINT*, "Pre-multiply the matrix of the score of the likelihood by its transpose"

 INFO = MATMUL(GHAT3PRM,GHAT3)


!----------------------------------------------------------------------------------------------------------------------
!(g) Invert the resulting num_param*num_param matrix
!    LINDS: compute the inverse of a real symmetric positive definite matrix
!    CALL LINDS (N, A, LDA, AINV, LDAINV)
!    N: Order of the matrix A. (Input)
!    A: N by N matrix containing the symmetric positive definite matrix to be inverted. (Input)
!       Only the upper triangle of A is referenced.
!    LDA: leading dimension of A exactly as specified in the dimension statement of the calling program. (Input)
!    AINV: N by N matrix containing the inverse of A. (Output)
!          If A is not needed, A and AINV can share the same storage locations.
!    LDAINV: Leading dimension of AINV exactly as specified in the dimension statement of the calling program. (Input)
!----------------------------------------------------------------------------------------------------------------------
 PRINT*, "Invert the resulting [num_param*num_param] matrix"

 CALL DLINDS(num_param, INFO, num_param, INVINFO, num_param)



!------------------------------------------------------------------------------------------------------------
!(h) Calculate standard errors and t-statistics
!    Degrees of freedom: N-theta_dim = 1426-90 = 1336
!    Two-sided test:                            80%     90%      95%     98%     99%    99.5%  99.8%  99.9% 
!                                               1.282   1.645    1.960   2.326   2.576  2.807  3.090  3.291      
!------------------------------------------------------------------------------------------------------------
 PRINT*, "Calculate standard errors and t-statistic"
 WRITE(3500,*) "Standard error of theta"
 
 m1 = 0

 DO m1 = 1, num_param
  
  SE(m1) = SQRT(INVINFO(m1,m1))
  WRITE(3500,*) SE(m1) 
   
 
 END DO

 WRITE(3500,*) ""
 WRITE(3500,*) "T-statistic of theta"
 

 DO m1 = 1, num_param
  
  TSTAT(m1) = theta(m1)/SE(m1)
  WRITE(3500,*) TSTAT(m1) 
 
 END DO


WRITE(3500,*) ""
PRINT*, "End of the computation of standard errors"
PAUSE  
 
!---------------------------------------------------------
END IF !if case for when the st errors are to be computed
!---------------------------------------------------------


DEALLOCATE(INVINFO, INVINFO2, INFO, INFO2)
DEALLOCATE(GHAT2, GHAT3)
DEALLOCATE(GHAT3PRM)
DEALLOCATE(SE, TSTAT,theta)

!****************************************
!(8) Likelihood maximization subroutine
!****************************************
WRITE(3000,*) "START THE AMOEBA COMPUTATION: "
WRITE(3000,*) ""
WRITE(3100,*) "START THE AMOEBA COMPUTATION: "
WRITE(3100,*) ""
PRINT*, "Start the AMOEBA computation: "
PRINT*, "here"



CALL amoeba(param_matrix,F,ftol,iter)

 DEALLOCATE(pos) !(ow CALL parameter_generation(X) within SUBROUTINE likelihood_function(X,F) cannot be executed)


!*********************
!(9) Output record
!*********************
OPEN(unit=777, file='track_like.txt',position="append")
!OPEN(unit=77, file='track_param.txt', position="append")
	
!WRITE(777,*) "Number of Function Evaluations ="
!WRITE(777,'(I11)') (iter)
WRITE(777,*) " "
WRITE(777,*) "Maximum Likelihood Vector ="
WRITE(777,'(F16.4)') (F)
!WRITE(77,*)  "Matrix of Parameters Maximizing the Likelihood:"
!WRITE(77,'(22F16.4)')   ((param_matrix(i1,12), i2=1,num_param), i1=1,num_param+1)

CLOSE(unit=777)
!CLOSE(unit=77)

CLOSE(3100)

!******************
!(10) Time record
!******************
!OPEN(unit = 2000, file = 'Time.txt', position = "append")
!WRITE(2000,*) (hour)
!CLOSE(unit=2000)
!CALL TIME(hour)
CALL cpu_time(etime)
WRITE(3000,*) "Program finished successfully (",(etime),")"
PRINT*, "Program finished successfully"
PRINT*, "Ending time", (etime)
!PRINT*, "Ending Time", hour

CLOSE(unit=3000)
CLOSE(unit=3500)

END PROGRAM main
